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We study flows and interface deformations produced by the scattering of a laser 
beam propagating through non absorbing turbid fluids. Light scattering produces 
a force density resulting from the transfer of linear momentum from the laser to 
the scatterers. The flow induced in the direction of the beam propagation, called 
"optical streaming", is also able to deform the interface separating the two liquid 
phases and to produce wide humps. The viscous flow taking place in these two 
liquid layers is solved analytically, in one of the two liquid layers with a stream 
function formulation, as well as numerically in both fluids using a Boundary 
Integral Element Method. Quantitative comparisons are shown between the nu- 
merical and analytical flow patterns. Moreover, we present predictive simulations 
dedicated to the effects of the geometry, of the scattering strength and of the 
viscosities, on both the flow pattern and the deformation of the interface. Theo- 
retical arguments are finally put forth to explain the robustness of the emergence 
of secondary flows in a two-layer fluid system. 
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1. Introduction 



Stre aming flow induced by light scattering has been emphasized by lSavchenko et al. 



(|l997T ) m various liquids such as liquid crystals and colloids. It was showed that 



a light beam having an inhomogeneous transverse intensity profile may enforce 
convection in liquids as a consequence of the transfer of linear and angular mo- 
mentum from light to liquids. 

In a recent inves tigation on the effe ct of a continuous laser wave on soft liquid- 
liquid interfaces (jSchroll et a/.ll2007l ). we proposed an original experiment to op- 
tically induce bulk flows in turbid liquids without any heating. We showed that 
light scattering by refractive index fluctuations can produce a steady flow due 
to transfer of momentum from light to the liquid. In addition, the viscous stress 
exerted by this flow can deform the liquid interface. We successfully compared 
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Figure 1. Eddies produced by the scattering of a continuous laser beam propagating in a non 
absorbing turbid liquid (micellar phase of microemulsion) layer bounded by a top glass wall 
and a bottom transparent brine solution; a thin layer of metastable foam is present at the 
brine-microemulsion interface. A beam centered toroidal eddy, evidenced by the advection of 
foam in this picture, induces a centered deformation of the interface by viscous withdrawal in 
the direction of beam propagation. The interface deformation is re presented by the bl ack curve. 
The experiment is performed in a so-called Winsor II equilibrium (Kella v et qZ.|[l994r j obtained 
here by mixing equal volum es of n-Heptane and brine with a small amount of AOT surfactant 
(|B. P. Binks fc P etsev 2000). Sodium chloride is used to screen electrostatic repulsion between 
surfactant heads and thus to significantly reduce the interfacial tension. For the chosen concen- 
trations, [NaCl] = 0.05M and [AOT] = 40mM, the interfacial tension is a = 1 - 5 10~ 6 N/m. 
The beam power is P — 1.32W and the beam waist luq = 2.99/xm. Eddies are magnified by 
overlaying 100 frames from a video sequence. 



experimentally observed deformations of the two-fluid interface to their analyt- 
ical prediction. We call this flow her eafter "optica l streaming" by analogy with 
the acoustic streaming phenomenon ( Nvborel 19581 ) that occurs when an acoustic 
wave propagates through a sound absorbing liquid. 

Indeed, while acoustic streaming is dissipative in nature due to wave absorp- 
tion, "optical streaming" results from the scattering of the incident beam on 
density inhomogeneities. However, the net result in terms of induced flow is very 
similar. More recently, we demonstrated that optical streaming is responsible 
for the droplet emission a t the tip of liquid jets triggered by radiation pressure 
(Wunenburg er et afJlioicj ). 

Optical streaming can be observed experimentally as illustrated in Figure [TJ 
In the present investigation, we consider theoretically the optical streaming tak- 
ing place in a set of two horizontal liquid layers separated by an interface. They 
are submitted to a laser beam propagating upward along the z-axis (see figure [2]). 
Assuming that at least one of the layers is a turbid liquid, i.e. contains refractive 
index inhomogeneities smaller that the optical wavelength, a bulk force is induced 
in this layer by light scattering leading to the formation of steady toroidal ed- 
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Figure 2. Sketch of the flow domain composed of the two fluids and the liquid interface 

initially at rest. 



dies. If the interfacial tension is sufficiently weak, these eddies are able to deform 
the liquid-liquid interface by the induced viscous stress. The goal of the present 
investigation is thus to study optical streaming in a more general frame, by an- 
alyzing the influence of the physical properties of the fluids and the geometry 
of the container. Beyond understanding the nature of these flows, an interestin g 
prospect is to use them in future microfluidic experiments ( Delville et al. 2009I ) . 
Our article is structured as follows: the physical modeling of the problem is pre- 
sented in section 2. In section 3, we analytically model the optical streaming in 
a single liquid layer, using a stream function formulation. The numerical simula- 
tion in a two layer system is presented in section 4, using a Boundary Element 
Method. In section 5, comparisons between numerical and analytical predictions 
are discussed, addressing the influence of the aspect ratio and sizes of the cell 
containing the fluids, of the scattering intensity and of the viscosity ratio on the 
flow patterns. Finally, good quantitative agreements between numerical simula- 
tions and experiments are reported in the appendix confirming the validity of the 
physical description and the numerical solution proposed in this work. 



2. Formulation of the problem 

2.1. Geometry and liquid properties 

We consider two liquid layers labeled by i, i=l ,2, separated by an interface, il- 
luminated by a continuous laser beam, of Gaussian intensity profile along the 
radial coordinate and of cylindrical symmetry along the vertical propagation 
axis z. The two liquids are enclosed in a cylindrical container of solid bound- 
aries S\yi and Sw2, of radius R ~ 10~ 3 m >> ujq ~ 10~ 5 m and thicknesses 
L2 — L\ ~ 10 -4 m >> ujq, where wn is the characteristic laser beam radius, also 
called beam waist (see figure [2]). Considering the axisymmetry of the system along 
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the z axis, we choose cylindrical coordinates (e r , ee, e z ) with origin O located at 
the intersection of the beam axis with the initially flat interface located therefore 
at z = 0. We refer to a point x by the space coordinates (r, 0, z). 
Physical properties of liquids i = 1, 2 are their viscosities r^, densities pi, refrac- 
tive indices ni and momentum attenuation coefficients Wi which represent the 
relative loss of light linear momentum per unit propagation length. The interfa- 
cial tension is denoted by a. The characteristic values of the liquids are rji ~ 10~ 3 
Pa s, pi ~10 3 kg m -3 , m ~ 1.33, W{ ~ 70m _1 and a ~ 10~ & N/m. 

2.2. Bulk force associated to light scattering 

The momentum flux associated to light of celerity c in vacuum, propagating in a 
fluid of refractive index rii is ^1. The beam intensity profile is 

2 P 

I(r, z) w 7(r) = ^2 e_2(r/tJo) (2- 1 ) 

P is the beam power. 

When light is elastically scattered by the refractive index inhomogeneities of the 
medium, the forward component of its momentum density qi z decreases along 
the direction of propagation. To quantify this decrease, we introduce the forward 
momentum attenuation coefficient in liquid i defined by 

ro . = _J_%* 2 (22) 

Qiz oz 

We assumed statistical isotropy of the distribution of the refractive index inhomo- 
geneities within the medium. Therefore, the total rate at which the momentum 
density decreases due to scattering is Wi^Ie z , i = 1,2. As a consequence of 
momentum conservation, the liquids are subjected to a body force density F,- L 
such as 

Ft = w^Ie z = F ie - 2 ^ 2 e z , i = l,2 (2.3) 
c 

2.3. Flow equations 

Since the characteristic order of magnitude of o ptically driven flows as m ea- 
sured in the experimental investigation reported in Wunenburger et all ( 20ld ) is 
Ui ~10~ 4 m s _1 , we investigate a viscous flow at low Reynolds number (Rei = 
UiLiPi/rn ~ 1(T 2 ). 

Both liquids satisfy the mass conservation and Stokes equations, 

V -u t = ; i = 1,2 (2.4) 

= -V Pi + TfcAtij + Fi ; i = 1, 2 (2.5) 

The gravitational force has been included in the corrected pressure term pi that 
is defined by pi = p[ + pigz, p\ being the fluid pressure. 

Assuming no slip at the interface and immiscible fluids, the hydrodynamic stress 
balance on the interface (described by its height h(r)) involving interfacial tension 
and gravity effects is written as 



B.C.I : T 1 -n-T 2 -n = (o-K-(p 1 -p 2 )gh)n on 5/ (2.6) 
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Here, Tj = —pil + ^(Vitj + Vu/) is the corrected hydrodynamic stress tensor, 
n is the unit vector normal to the interface directed from fluid 1 to fluid 2 and 
1 d r— 

n(r) = = is the double mean curvature of the axisymmetric interface 

in cylindrical coordinates. The continuity of velocity at the interface yields 

B.C. 2 : u 1 = u 2 on 5/ (2.7) 

At the container walls, the classical no-slip condition for viscous flow is assumed, 
i.e. 

B.C.3 : U! = on S W i \ B.CA : u 2 = on S W 2 (2.8) 
Finally, the motion of the interface is described using a Lagrangian approach 

dx 

— = u(x) on 5/ (2.9) 



3. Analytical solution of a one-fluid model 

When the top fluid is much less viscous and much less dense than the bot- 
tom one, a simplified analytical solution to the above problem can be derived. 
Indeed, in this situation of an inviscid top fluid remaining at a constant uniform 
pressure, the stress exerted on the bottom layer can be neglected leading to a 
one-fluid model. To develop the solution, the flow taking place in the bottom 
fluid layer can be solved analytically using a stream function formulation. 
We assume that the width of the laser beam is much smaller than the width of 
the resulting eddy (in figure CD for instance, the eddies presented are ten times 
larger than the beam waist). Therefore, the scattering force can be represented 
as a line force f\e z = F\-Kuj\e, z along the centre line of the cylindrical container. 
The amplitude of the deformation is supposed to be very small compared to its 
characteristic width and a flat interface is assumed. 



The flow within the layer is modeled as a steady viscous flow contained in a 
cylindrical layer of height L\ and radius R centered on the beam axis. Due to 
the symmetry of the optical excitation, the flow is axisymmetric, and can be 
described by a stream function ^>(r, z) such that 

u = V x -^^e e (3.1 

r 



This choice for the stream function ensures incompressibility (Equation ([27 
and the Stokes equation becomes 

= -Vp + ??iV 2 (V x t^l e@ ^ +Fl (3.2) 

Taking the curl of Equation (|3.2p removes the pressure gradient, and yields 

n 1 D A ^(r,z)=r(VxF 1 )-e @ (3.3) 
where D 4> = D 2 (D 2 -) and D 2 = J^j — ^-§p + Assuming a line force, the 
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right-hand term of Equation (|3.3p is a delta-function at r = 0. Thus, except at 
r = 

L>V(r,z) = (3.4) 
Assuming separation of variables, the general solution is 

ip (r, z) = G(r)sin kz (3-5) 

where k = n/L\. 

Note that this solution always satisfies ^f(r, z = 0) = and ^-(r, z = — Li) = 0. 
Even though not strictly correct, this last condition ensures a straightforward 
analytical solution. 

The solution to the problem is finally achieved using four boundary conditions. 
We require by symmetry u r (0, z) = 0, no flux through the far wall of the container 
r = R, u r (R, z) = and no slip along it u z (R, z) = 0. 
Replacing Equation (|3.5p into Equation (j3.4|) we find 

jj / fkr fkr 

G (r) = —r I K\ (kr) / I\ (x) K\ (x) xdx — I\ (kr) \ K\ (x) 2 xdx 

4vrfc V JkR JkR< 

tj / fkr fkr 

f rA(k) ( h (kr) / h (x) K x (x) xdx - K x (kr) \ h (xfxdx ](3.6) 

V JkR JkR 



Airk 



where A(k) = ^° Jr^^I^^ \ h an d Ki being the first order modified Bessel 

J () Ii(x) xdx 

functions of first and second kind respectively. 

The last free parameter Ut can now be estimated by matching the flow amplitude 
and the forcing. This means that the integrated delta-function forcing must be re- 
lated to the scattering induced flow. More specifically, we take the cross-sectional 
integral of Equation (|3.3p at the mid-plane of the liquid layer, i.e. at z = —L\/2. 
On the left-hand side, we get 

771 J D 4 iidS = rjiUf (3.7) 
On the right-hand side, combining Equations (|2.ip and (I2.3p . we get 
J r (VxFi) • e e dS = 2n ' r 2 (-^T\ dr = 2Pvo x ^ = h 



(3.8) 



Thus, the characteristic velocity of the light-induced bulk flow Uf is 

Uf=^=^ (3.9) 



The velocity components are deduced from Equations (|3.ip and (|3.5p 

G(r) 

u r (r,z) = kcoskz (3.10) 



fkr fkr 

-Ko(kr) / I\ (x) K\ (x) xdx — Iq (kr) / K\(x) 2 xdx 

JkR JkR 



sin kz 
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jj / rkr rkr \ 

— — —A(k) ( Iq (kr) / I\ (x) K\ (x) xdx + Kq (kr) / I\ (x) 2 xdx ) sin kz 
4?r V JkR JkR J 

(3.11) 

Assuming p(r, z = —Li/2) = 0, the corrected pressure field can be deduced from 
equation (3^ by noticing that J = -r/i±^ and §f = r?^^- It yields 

p(r, z) = rjt-^k [K (kr) + A(k)I (kr)] cos kz (3.12) 

where Jo an d Kq are the 0* ft order modified Bessel functions of first and second 
kind respectively. The solution of this simplified model allows to express the 
corrected normal stress at the top of the fluid layer which is given by the T zz 
component of the corrected stress tensor 

T zz (r,z) = - P (r,z) + 2 m -^ (3.13) 

Using Equations (13. lip and A3. 12H . we can write 

Ut 

T zz (r,z) = —7]i—kS(r)coskz (3-14) 

7T 



with 



S (r) = v " ' (1-2 r h (x) K x (x) xdx) - f Kl (xfdx 



kR J L JkR 



A(k) 



(pkr \ pkr 

1 + 2 h(x) K x (x) xdx I + 2K (kr) / I x (xfxdx 
JkR, J JkR. 



(3.15) 

To predict the deformation, the corrected normal stress exerted by the viscous 
flow in the lower layer is balanced with buoyancy and capillarity. Since the in- 
terface deformation is considered to remain small (i.e. h! = dh/dr « 1), the 
corrected hydrodynamic stress experienced by the interface at radius r is well ap- 
proximated by —T zz (r,0)e z . The linearized equilibrium equation of the interface 
is therefore written in a dimensionless form using L\ as a characteristic length 
and <jjL\ as a characteristic pressure 

h"(f) + ^ - Bo h(r) + ^S(f) = (3.16) 
r a 

where Bo = p\gh\ja is the gravitational Bond number and the symbol 7 refers 
to a dimensionless quantity. 

This equation is solved using a differential algebraic equation solver (DASSL) 
based on a backward differentiation formula in the interval [r = e,r = R]. Note 
that e ^ since T zz (0, 0) diverges due to the line force distribution exerted at r = 
0. Since Equation (|3.16|) is of second order, we require two boundary conditions 
for h(r). We therefore assume h(r = e) = h t and need the condition dh/dr(r = 
e) = h' a h e being our shooting parameter. The convergence condition is chosen 
to be dh/dr(r = R) = with a tolerance of 0.1%. Analytical profiles predictions 
are compared with experimental and numerical results in the appendix. 
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To check the predictibility of the present simplified one-layer analytical model, 
we also investigate numerically the two-fluid configuration. This is presented in 
the following section. 



4. Numerical resolution 



The Boundary Element Method reveals to be very accurate for solving inter- 
facial flow problems with high resol ution as reported in the analysis of flow in- 
volving electric and magnet ic fields ( Sherwoodlll987h. optical rad i ation pressure 
dChraibi et al\ 1200861 . l2010l ) or buovancv (|Manga k Stonel ll994J : iKoch fc Koch 
19941 ). When used to solve problems involving surface forces like capillarity, radi- 
ation pressure, etc.. and volume forces, which can be written under a conservative 
form, this method allows to reduce the formulation of the problem to boundary in- 
tegrals only. In the present case, the scattering force obviously does not belong to 
this class, and consequently the integral formulation of the proble m involves a vol 



ume i ntegral term correspond ing to the contribution of this force (|Occhialini et al 



1992!; Ilssenmann et ad 1201 ll ). The integral formulation is achieved by first mak- 
ing use of the Green's function for the Stoke s operator which corresponds to the 
solution of the fundamental Stokes problem ( Pozrikidisl 1992 ). The integral equa- 
tions are presented below under a dimensionless form by choosing U a = cr/r]\ as 
a reference velocity, L\ as a reference length and pio = rua 7 '{L\r)\) i = 1,2 as 
reference pressures in both fluids. Taking into account the boundary conditions 
on the interface Si and on the walls Swii Sw2 (see figure [2j), and denoting by 
the outward unit normal vector to these boundaries, the integral form provid- 



n 



ing the velocity at any point x in the interior of the volumes V\ and V 2 , is given by 



u 



:(*) = [ 
Js 



U-Ti-n e dS v 



Si+Swi 



/ n e -K-UidS y + / U-FidVy i = 1,2 

J Si+Swi JVi 

(4.1) 

The solution of this equation requires a prior determination of stress on Sj, Swi 
and Sw2 as well as of the velocity on Si. This is achieved by solving an integral 
equation which, once B.C.I (equation (|2.6p ) is explicitly taken into account, can 
be written as 



i + C 



-u[x 



/ U • n (h(fy) - Bo h(f y )) dS y + 

JS! 



(C - 1) / n • K • u dS y + / U-Tx-ndSy-C U -T 2 -n dS y + 

J Sj J Swi J Sw2 



U-F 1 dVy + U U'F 2 dVy 

Vi Jv 2 



(4.2) 



where x locates a point on Si, n = riyi and C = V2/V1 is the viscosity ratio. 
Her e U and K are G reen kernels for velocity and stress respectively and are given 
by rtPozrikididli992r ) 



(4.3) 
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, , 3 ,ddd, 
K{d) = -^ { ^ * = 1 ' 2 (44) 
where d = y — x, y(r y ,z y ) is the integration point. The first term in the right 
hand side of equation (14. 2p describes the flow contribution from interfacial ten- 
sion and gravity, whereas the second term accounts for shear rates contrast on the 
interface. This term vanishes when £ = 1. The third and fourth terms account for 
shear occurring on Syyi an d Sw2 as a result of the no-slip boundary condition. 
The two last terms account for the contribution of the flow due to the scattering 
force density. Unlike the other terms, they are evaluated on the volumes V\ and 
V2 respectively for fluid 1 and 2. Velocities on the interface as well as stress over 
all the boundaries Si, Swi and Sw2 are determined by solving the discrete form 
of this equation using a numerical pro cedure. This is performed by an analytical 
integration in the azimuthal direction ( Lee &: Leal 19821 ; Graziani 19891 ) reducing 
all surface boundaries Si, Swi, Sw2 into line boundaries Ti, Tw\, ^W2 and vol- 
umes V\ and V2 into surfaces Si and S2 (see figure [2]) . The numerical procedure 
requires then the discretization of all line boundaries and surfaces with line seg- 
ments and quadrilaterals respectively. Line discretization makes use of constant 
boundary elements, i.e. line segments with centered nodes. Elliptic line integrals 
resulting from the azimutal i ntegration are evaluated using Gauss quadratures 
and power series expansions (|Bakrlll985l ). Numerical integr ation on quadrilater- 
als is performed using two-dimensional Gauss quadratures ( Davis h Rabinowitd 
1984 ) and isoparametric transformations. The fluid-fluid interface is parameter- 



ized in terms of arc length and is approximated by local cubic splines, so that 
the curvature can be accurately computed. Distribution and number of points 
are adapted to the shape of the interface, so that the concentration of elements 
is higher in regions where the variation of curvature of the interface is larger. 
The motion of the interface is followed using the kinematic condition (I2.9P which 
is discretized using an explicit first-order Euler time scheme. A typical computa- 
tion begins with a flat interface at rest. The laser beam is switched on at t=0, and 
the deformation of the interface begins. Computation stops when an equilibrium 

state is reached (— — > 0). The time step is chosen to be about 50 times smaller 

than the reference time r = L\/U a . 



5. Results 

5.1. Single layer flow 

In this section, we compare the steady velocity field in the bottom layer predicted 
using both the numerical simulation and the analytical one-layer model. For this 
comparison, we chose 772/^1 = 10~ 3 and p2 << pi while /2 = 0. Moreover, the 
influences of the layer aspect ratio R/L\ and of the layer thickness to beam size 
ratio Li/ujq on the flow field and interface shape are discussed. 

5.1.1. Flow pattern 

Comparisons between the flow patterns predicted analytically and numerically 
are reported in figure [3j It shows the r and z dependences of the axial (figures 
Wib-d)) and radial velocity components (figures ^ a, e-g)) in the bottom layer. 




Figure 3. Comparison between numerical (symbols) and analytical (lines) velocity fields in the 
bottom layer. The left figures show the variation of the reduced velocity components u r /Uf and 
u z /Uf with respect to r/L\ (Uf = /i/^i). The right ones present the variation of u r /Uf and 
u z /Uf with respect to z/Li. fx/a = 0.157, fi = 0, ^2/7/1 = 10 -3 , R/Li = 5, L\ = L2 and 
Li/ojo = 10. In the numerical resolution Bo = 5 and therefore h(r = 0)/L\ = 0.008. 

A satisfactory agreement is observed concerning u z (figures \3[b-d)) whatever the 
value of z. 

The same remark holds for u r (figures [3]( a,e-g)) at the horizontal mid-plane of 
the layer (z = —L\/2) and at the fluid-fluid interface (z = 0). A discrepancy 
can be identified at the bottom wall (z = —Li) due to the no-slip condition 
(u r (z = —L\) = 0) in the simulation while a free slip condition was chosen in 
the one-fluid model for simplicity purpose. For all radial velocities we can see 
that u r > for z > —L\/2 and u r < for z < — Li/2, which is consistent with 
the statement of a toroidal eddy centered vertically on the mid-plane of the fluid 
layer. 

Even though the one-layer model is a simplified version of the real situation, the 
comparison shows that this analytical model captures the essential trends of the 
flow. In the next section, we investigate the effects of the aspect ratio R/L\. 
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Figure 4. Axial velocity for several values of RjL\. fi/a = 0.157, fo = 0, 772/171 = 10~ 3 , 
L\ — 1/2 and Li/uq — 10. In the numerical simulation Bo — 5. Open symbols represent the 
numerical simulation while continuous lines represent the analytical resolution. Left inset shows 
the variation of u z (r = 0.05)/t// versus R/L\. Right inset shows the variation of r(u z = 0)/Li 
versus R/L%. 



5.1.2. Influence of the container aspect ratio 

Comparisons between the steady axial velocity profiles u z , predicted analyti- 
cally and numerically, for different values of the aspect ratio R/Li, are reported 
in figured! It shows a good agreement between the simulations and the analytical 
model except for |r|/Li < 0.1 i.e. within the beam |r| < ujq, since L\/ujq = 10. 
This result was expected since the analytical solution shows a divergence of u z 
for r = as a consequence of the line force distribution approximation. 
The left inset of figure H] shows that u z scales like Ut \og(R/ L\) for narrow tubes 
(L\ » R) while it scales like Ut in the thin films case (R >> L\). The right in- 
set of figure H] shows that the radial extension of the flow r(u z = 0) scales like the 
smallest dimension of the container. In fact, when L\ >> R, the flow is bounded 
by the lateral walls and consequently its radial extension scales like R. In the 
other limit, when R >> Li, the flow is unbounded in the r direction and mass 
conservation in the axisymmetric geometry produces toroidal eddies of circular 
section of radius L\. Deformation profiles predicted by the numerical simulations 
for different values of R/L\ are presented in figureEJ The left inset shows that the 
tip deformation scales like L\ when R >> L\. Similarly to the velocity field, the 
radial extension of the deformation scales like the smallest container dimension 
as shown by the right inset of figure [5l To conclude, in the limit of narrow beams 
[L\ >> wo), the radial extensions of the flow and the deformation always scales 
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Figure 5. Numerical deformation profiles for several values of R/L\. fi/a = 0.157, /2 = 0, 
772/771 = 10~ 3 , Bo = 5, Li = 1/2 and Li/cjo = 10. Left inset shows the variation of h(r = 0)/Li 
versus R/Li. Right inset shows the variation of r(h = 0)/Xi versus R/Li. 

like the smallest dimension of the layer. In the next section, the influence of the 
beam size on the flow and on the interface deformation is investigated. 

5.1.3. Influence of the beam size 

Comparisons between the steady axial velocity profiles u z , predicted analyti- 
cally and numerically, for different values of the layer thickness to beam waist 
ratio Li/ojq, are reported in figure [6j It shall be noted that the analytical profile 
is unique since the solution does not depend on ojq. We observe that all numerical 
profiles coincide for |r| > ojq. In addition, we notice a disagreement for \r\ < Wo, 
similarly to the results presented in figure HI While the analytical prediction shows 
a divergence of u z as r — > 0, the numerical simulation shows a saturation. This 
saturation always occurs at r ~ ujq. We conclude, in the limit of thin films, that 
there are two characteristic dimensions for the axial velocity. When |r| < o;o, 
the axial velocity depends on coq while its radial extension depends on the layer 
thickness L\ (see right inset of figured]). The left inset of figure [6] shows that u z 
scales like Uf iog(Li/wo) in the limit of thin films and narrow beams. A similar 
behaviour has been observed, in the Re « 1 limit, in a theoretical investiga- 
tion of the acoustic streaming velocity field by Nyborg (Hamilton fc Blackstock 
19981 ). Indeed, it was demonstrated that the maximum axial velocity of acoustic 
streaming increases logarithmically with the focal length (equivalent to L\) of 
the acoustic beam propagating through an infinite fluid. 

Deformation profiles predicted by the numerical simulations for different values 
of Li/ujq are shown in figure [3 We conclude from the insets that both the defor- 
mation height and its radial extension scale like L\ when L\ >> ujq. The general 
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Figure 6. Axial velocity for several values of Li/wo- f\/o = 1-57, fa = 0, 772/771 = 10 3 , L\ = L2 
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Figure 7. Numerical deformation profiles for several values of Li/ujq. fi/cr — 1.57, fa = 0, 
V2/V1 = 10 -3 , -Bo = 5, Li = L2 and R/Li = 5. Left inset shows the variation of the tip 
height h(r — 0)/Li versus L\/ujq. Right inset shows the variation of the radial extension of the 
deformation r(h — 0)/L\ versus L\/ujq. 
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conclusion of this section is that for thin films R » L\ and of narrow beams 
L\ >> a;o, both the radial extension of the axial velocity and of the deformation 
scale like the layer thickness L\ while their amplitudes scale respectively like Uf 
and L\. These limits correspond to our experimental investigations reported pre- 
viously (|Schroll et all 120071 : IWunenburger et fllll20inh . 

The following sections present the hydrodynamic interaction between the two 
layers with an investigation of the influence of the different aspect ratios R/L\, 
L1/L2 and force and viscosity relative amplitude /1//2 and 772/771- 
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Figure 8. Steady flow pattern for R/Li = 5 and Li/wo = L^/ojq — 30, S = /1//2 = 1, 
C — nijr\\ = 1, Bo — 10 and = 0.3. Velocity magnitude is reduced by U a = cr/771. 



5.2. Double-layer flow 

5.2.1. i<7ow patterns 

The flow patterns in wide containers (R > L\,L\ = L2) and long containers 
(R < L\ : L\ = L2) are displayed in figure [8] and figure El^a). 
We first observe in figure [HJ nearly circular streamlines, which are the signature 
of toroidal flows also called eddies. These eddies exist at steady state because 
(i) the flow has to follow tangentially the interface due to the continuous axial 
forcing of the scattering force and (ii) mass conservation induces a circulation 
of the fluid. When R > L\, figure [8] clearly shows that the spatial extension of 
these eddies scales like L\ as discussed in the previous sections (figure H|). When 
L\ > R (figure a)), the eddies invest the whole thickness of the layers while 
their width equals the container radius. To conclude, these results demonstrate 
that the axial dimension of the eddies is always L\ while the radial one scales 
like the smallest dimension of the fluid layer. 

Figures EJfr-c) display the flow pattern when L\ and L2 are different. We notice 
first that the maximum velocities are different in the layers despite an equal scat- 
tering force density. This is due to the dependence of the velocity with respect 
to the thickness of the fluid layer, which has been evidenced in the left inset of 
figure [H It should be also noticed that in the fluid layer of aspect ratio L\ / L2 
equal to 1, there is an eddy extending to the whole layer, while in the other layer 
of smaller height, a secondary eddy develops. In fact, because the eddy is larger 
in the fluid layer of largest thickness, it forces a motion in the other layer, as 
a result of momentum transfer by viscous shear through the interface, yielding 
to a secondary eddy. The scattering force exerted in both layers inducing two 
contrarotative eddies, any loss of symmetry of the flow domain induces the ap- 
pearance of a secondary corotative eddy within one of the layers. Similar patterns 
can be obtained by imposing non equal scattering force densities in the two layers 
as illustrated in the next section. 



16 H. Chraibi, R. Wunenburger, D. Lasseux, J. Petit and J-P. Delville 




Figure 9. Steady flow patterns, (a) R/Li = 2/3 and Li/cjo = L 2 /lo = 30 (6) R/L\ = 3, 
Li/ojo = 10 and L2/0J0 = 30 (c) R/Li = 1, Li/ljq = 30 and = 10. In all cases /i/a = 0.3, 

E = /1//2 = 1, C = Vz/ni — 1 an d Bo = 10. Velocity magnitude is reduced by t/ CT . 



5.2.2. Influence of the scattering force ratio 

In the present section, we discuss the dependence of the flow and of the inter- 
face deformation on the magnitude of the scattering force. 

Figure [10] shows the variation of h max = h{r = 0) with the total scattering 
force / = /1 + /2- We first notice that when the forces are equal in both 
layers fx = fz = //2, the variation of h max is linear with respect to /. For 
[fi = 0, /2 = /] and [fi = f, /2 = 0], we observe a small deviation from the linear 
behaviour, which is attributed to the dependence of h max to the layer thickness as 
shown in figure [71 Considering the case [/1 = 0, f'2 = f], the layer thickness above 
the hump (Lq — /12) is a little smaller than Lq leading to a smaller deformation, 
while in the case [f\ = /, /2 = 0], the layer thickness below the hump (Lq + h\) is 
a little larger than Lq leading to a larger deformation (see left inset of figure [7J . 
The right inset of figure [TO] shows that the variation of the difference of heights 
Ah = h\ — h>2 between the case [/1 = /, /2 = 0] and the case [/1 = 0, /2 = /] is 
quadratic with respect to /. 

In the left inset of figure [TOl we can see that the maximum velocities vary 
linearly with /, and that the same velocity is obtained in both configurations 
{[fx =f,f2 = 0] and [fi = 0,/ 2 = /]). In the configuration h = h = f/2, the 
maximum velocity in each layer is half that obtained in the other configurations. 
Consequently, the effects of the flows on interface deformations are additive as 
each layer induces its own viscous stress on the interface. 

Considering the flow pattern, we observed in figure [8] that when E = /1//2 = 1, 
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Figure 10. Hump height variation against the total scattering force / when [fi = f,f 2 = 0], 
[/i = 0, h = f] and /j = f 2 = //2 (7?/Li = 3.33, Li/u = L 2 /u;o = 30, Bo = 45 and 771 = 772). 
The dashed line is a linear regression for fi = f 2 = f/2. The left inset presents the maximum 
fluid velocity in each configuration while the right inset shows the variation of the reduced 
difference of heights Ahmax/Li between the cases [fi = f,f 2 = 0] and [/1 = 0, f 2 = f] against 
f/o. 



two co-rotating eddies are induced, both of them rotating clockwise, as a result 
of the direction of the scattering force density along the z axis. In this perfectly 
symmetric case, no secondary eddy is generated. 

In figures [TTJ [12] and [131 corresponding to E / 1, a secondary eddy rotating 
counter-clockwise develops in the fluid layer where the amplitude of the scatter- 
ing force is the smallest. This secondary eddy, whose extension increases with E, 
is contra-rotative compared to the main eddies. 

Finally, in figure [U corresponding to [fi = /, /2 = 0], we observe two contra- 
rotating eddies, the one in the bottom layer being determined by the direction of 
propagation of the laser beam while the one in the top layer, which is of weaker 
intensity, is exclusively induced by viscous shear through the interface as a result 
of zero scattering force density, = 0. 

To summarize, when / < a, there is no difference whether the scattering force is 
applied in the top or in the bottom layer. When / > a, a small difference between 
the deformations can be noticed and is attributed to the resulting difference of 
thickness of the scattering layer. The flow patterns showed that a difference of 
scattering force between the layers leads to the emergence of a secondary eddy 
in the layer where the force is the smallest. 

Up to now, we have investigated flow patterns for £ = 772/771 = 1 in order to 
analyze the influence of the scattering force amplitude and of the different aspect 
ratios of the container. However, since viscosity affects flow velocities, the ques- 
tion remains on whether it modifies the two-layer flow pattern driven by light 
scattering. This is the purpose of the following section. 
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Figure 11. Steady flow pattern for R/Li = 5 and Li/uo = L 2 /ujo — 30, £ = /1//2 = 2, 
C = 7)2/771 = 1, Bo = 10 and fi/a — 0.3. Velocity magnitude is reduced by U a . 




Figure 12. Steady flow pattern for R/L\ — 5 and Li/cj = £2/1^0 — 30, S = /1//2 = 4, 
£ = 772/771 = 1, Bo = 10 and fi/u ~ 0.3. Velocity magnitude is reduced by [7 CT . 



5.2.3. Influence of the viscosity ratio 

In figure [15l the variations of the maximum velocity u max , in fluid layer i is 
represented versus £ = r/2/771 while keeping S = /1//2 = 1 and Li = Li- 
As we set 7/1 to a constant value, there is no variation of the velocity magnitude 
in layer 1, while the velocity in layer 2 scales like l/£ = T7i/f?2- This demonstrates 
that the dimensional scaling law established for a single layer flow U{ oc FiUj^jrn 
(fi = Fi-KoJ^) is valid in each fluid layer as wo is a characteristic length for the 
axial velocity near the z-axis (see comments of figured]). Moreover, the inset of 
figure US] shows that h max does not depend on the viscosity ratio. This can be 
explained by recalling that T zz i is proportional to r/jitj while Ui is proportional to 
1/rji making T zz i and hence h max independent of the fluid viscosities. 
In order to get a general understanding of the location of the secondary eddy, a 
parametric study of the flow pattern was performed by varying £ and C, in the set 
(S, C) S {0.01; 1; 100} 2 while keeping the total scattering force / = fx + fa con- 
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Figure 13. Steady flow pattern for R/Li = 5 and Li/wo = £2/^0 = 30, E = /1//2 = 1/4, 
£ = 7/2/771 = 1, -Bo = 10 and fi/a = 0.3. Velocity magnitude is reduced by C/ CT . 




Figure 14. Steady flow pattern for R/L\ = 5 and Li/wo = £2/^0 = 30, /i/cr = 0.3, /2 = 0, 
C = n'ljnx — 1, Bo = 10. Velocity magnitude is reduced by U a - 
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Table 1. Secondary eddy location for (E,C) 6 {0-1; 1; 10} 2 . Bottom layer, top layer or 
inexistent secondary eddy are respectively marked as 1, 2 and 0. 
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stant. Table Q] summarises these results, indicating whether the secondary eddy 
is in the bottom layer (1), top layer (2) or nonexistent (0). 

Table Q] shows that the location of the secondary eddy is independent of Q and 
is always located in the layer were the scattering force density is the smallest. 
This can be explained by noting that, since F% oc tjiUi/ojq and F2 cx ^vq/oJq, one 
has u\/v,2 oc ££. This scaling law is confirmed by our computational results, as 
illustrated by the two following cases. 

In the first case, X = 10 and £ = 0.1, we found that the maximum velocities in 
layers 1 and 2 are u max 1 

— u max 2/Ufj — 0.012. In the second case, S — 0.1 
and C = 0.1, we found u max 1 /U a ~ 0.00114 and 

^max 2/U a — 0.114, in agree- 
ment with the prediction u\/u2 = 0.01. The viscosity ratio only determines the 
velocity magnitude in each fluid layer. This observation is confirmed by a simple 
model of parallel flow which is presented in the following section. 



5.2.4. Lubrication model 

The location of the secondary eddy is predicted by making use of a lubrication 
model assuming R >> Li >> ujq and derived from a two-dimensional cartesian 
formulation of the Stokes equation. For the sake of simplicity, we consider a two- 
dimensional flow of velocity u = (u x , 0, u z ) in cartesian coordinates (x, y, z). Since 
we are only interested in the distribution of eddies i.e. in the z dependence of u x , 
we consider layers of infinite horizontal extent and a parallel flow u = u x {z)e x . 
Therefore far enough from the z-axis along which the forcing is exerted (x >> 0), 
we can write the Stokes equation in both layers, as follows 

(«)0 = -^ + TH-MT ; (6)0 = -— ; , = 1,2 (5.1) 
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We can deduce from (|5.1I (b)) that the flow is driven by a horizontal pressure 
gradient. 

Differentiating equation (|5.1I (a)) with respect to z, we find that 

= ^H; i = l,2 (5.2) 

This leads to the following expression of the velocity in each layer 

u ix (z) = a,iz 2 + biZ + Ci ; i = 1, 2 (5.3) 

where a^, 6j, q, i = 1,2 are six numerical constants to be determined by using the 
six following conditions 

ui x (z = 0) = u 2x (z = 0) (5.4) 

duix I n\ du2x I n\ (k k\ 

m — (x = 0)=m- w (z = 0) (5.5) 

o 

u\ x (z)dz = (5-6) 



L2 







u 2x (z)dz = (5-7) 
-Ft (5.8) 



(5.9) 



dx 

— = F 2 

dx 

The first two conditions (|5.4p and (|5.5p simply express the continuity of velocities 
and shear stress along the interface. Conditions (|5.6p and (|5.7p express mass 
conservation in each layer. Conditions (|5.8p and (|5.9p express the fact that the 
pressure gradient driving the flow is related to the scattering force driving the 
motion near the explained hereafter. Because the eddies always scale 

like L{ when R » Li, we consider that ^ oc 3^ cx Pi/Li. In addition, the 

Stokes equation near the z-axis leads to ^ oc rji oc Fi and finally to ^ oc 
Pi/Lj oc Fi. One can also notice a (-) sign on equation (|5.8p . This is justified by 
the fact that when the system is symmetric (i.e. Fi=F 2 , rji = rj 2 , L\ = L 2 ), we 
assume that there is no motion of the interface along the x axis, which leads to 
opposite pressure gradients and velocities near the interface {u\ x (z = CP) > 
and u 2x (z = + ) < 0). The solution approximated the flow without satisfying the 
no-slip boundary conditions at the top and bottom walls. The conditions imposed 
to the two-fluid, parallel flow lead to 
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Figure 16. Velocity u x (z) predicted by the cartesian model, (a) Symmetric case: E = 1, £ = 1 
and L1/L2 ~ 1. (b) Contrast of magnitude of scattering force densities: S = 2, £ = 1 and 
L1/L2 = 1. (c) Contrast of viscosities: E = 1, ( = 0.2 and L1/L2 = 1. (d) Contrast of layer 
thicknesses: E = 1, f = 1, -L1/L2 = 1.5. 



«l = 5.12 

a 2 = -i (5.13) 
2r/ 2 

Figure [16] illustrates four interesting cases highlighting the emergence or not of 
a secondary eddy in one of the layers. In figure [T6T a), the system is symmetric 
and u x (z) changes sign only once in each layer. This means that only one eddy 
exists in each layer. In case (b), the magnitude of the scattering force density in 
the top layer is half that in the bottom one. In the top layer, u x (z) changes sign 
twice (at z = 0.79 and z = 0.21). This means that two contra-rotative eddies 
exist in the top layer. In case (c), the fluid in the bottom layer is 5 times more 
viscous than the one in the top layer. We can notice that although the velocities 
in the bottom layer are much smaller, there is no secondary eddy. Finally, in 
case (d) the bottom layer is 1.5 times thicker than the top layer. A secondary 
flow appears in the layer of smallest thickness, since u x (z) changes sign twice at 
z = 0.73 and z = 0.14. As a conclusion, even though very simple, this parallel 
flow model gives a qualitative prediction of the location of the secondary eddy 
confirming the numerical simulations of the axisymmetric flow. 
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6. General discussion 

Optical streaming represents a new way to drive fluid flows by light. It can be 
useful in situations where the direction and magnitude of fluid flows need to be 
controlled. The only condition for inducing these optical streaming flows is to 
use non absorbing turbid fluids. Any viscous fluid can be considered, while the 
strength of the optical streaming will strongly depend on the s ize and nature of the 



scatte r ers which could be s olid or liquid nanosuspensions (see IWunenburger et al. 



temnMSchroll et all (|oo3) "or an example of density fluctuations in near-critical 



binary mixtures acting as liquid scatterers) and generally speaking any fluid pre- 
senting refractive index inhomogeneities at a length scale much smaller than the 
optical wavelength. 

Microfluidics is of particular interest for application of optical streaming since the 
small dimensions involved are compatible with the laser focusing. If a fluid-fluid 
interface is present in these experiments (either microfluidics coflow experiments 
or digital microfluidics), we expect the emergence of secondary eddies depending 
on the thicknesses of the fluid layers and on the scatterers present in these flu- 
ids. In addition, if the interfacial tension is sufficiently low, this could also lead 
to interface deformation. While a viscosity contrast will not influence the loca- 
tion of secondary flows, a significant impact of this contrast on flow magnitude 
is expected since velocity is inversely proportional to viscosity. Many applica- 
tions of optical streaming can thus be devised. For instance, flow control at very 
small flow rates is always difficult in microchannels as syringe pumps operate 
close to their mechanical limits; flow rate changes are also always accompanied 
by complex transients. In these regimes, the fine and smooth control of viscous 
flows by optically tunable scattering forces is probably an interesting alternative 
for turbid fluids. Mic romixing is also very challenging in classical microchannels 
(INguyen because flows are mostly laminar and thus constitute a 



major drawback for efficient chemical reactions (jBaroud et adl2003l ). One route 
for acti ve mixing consi sts in superimposing pressure perturbations to the main 
stream ( Lee et al. 200 ll ). Scattering forces could replace pressure actuation if flu- 



ids are turbid or seeded with appropriate scatterers. 



7. Conclusion 

In this work, we investigated flow production by light scattering in a two-fluid 
system involving turbid fluids. We compared the flow pattern predicted numeri- 
cally to the analytical one-fluid model prediction. This simplified analytical de- 
scription predicts most of the characteristic features of both the induced flow and 
the resulting interface deformation. Due to the large scale of induced flows, we 
investigated the geometrical effects of the container. We showed that the veloc- 
ity field and the interface deformation scale like the smallest dimension of the 
container. Moreover a variety of flow patterns were illustrated when breaking the 
system symmetries. For instance, in absence of scattering force intensity contrast, 
secondary eddies always emerge in the thinnest layer while in absence of layer 
thickness contrast, secondary eddies emerge in the fluid layer where the scattering 



24 H. Chraibi, R. Wunenburger, D. Lasseux, J. Petit and J-P. Delville 

force density is the weakest. On the contrary the variation of the viscosity ratio 
only affects velocities within the fluids while neither the interface deformation 
nor the flow pattern at steady state are influenced. Finally, a lubrication model 
was used to explain the emergence and to predict the location of these secondary 
flows in the layers. 

Generally speaking, light-induced flows exist whenever isotropic scattering flu- 
ids are at work; fluids near a second-order phase transition can be an interest- 
ing example to perform experiments. Consequently, we have proposed and well- 
identified a general non dissipative way to produce bulk flow with light that goes 
one step forward in the exploration of the coupling between light and hydrody- 
namics, so called " optohydrodynamics" . 



8. Appendix : Comparison with experiments : Interface deformation 
in the presence of radiation pressure 

In this appendix, we show comparisons between numerical, analytical and ex- 
perimental steady profiles of interface defor mations. In experim ents involving 
turbid two-phase systems like the one used in ISchroll et al\ (j2007l ) , the interface 
is made visible thanks to a refractive index contrast. In addition to the hydro- 
dynamics stress due to the op t ical s treaming, radiation pressure is exerted on 



the interface (see IChraibi et all l|20ld )). In this append ix, both stres s es are taken 



into account in the numerical simulation. Note that in lSchroll et al\ (120071). onl y 



(j20ld ). 



optical streaming effects were taken into account while in lChraibi et al 
only radiation pressure effects were modeled. 

In order to model radiation pr essure effects, we ha ve to add the following term 
(which is presented in details in lChraibi et al. ( 20ld )) in the right side of equation 



n(r) 



I(r) 



n 2 cos 0;(2(1 - 5) cos (9, - T(0i,0 t )((l - 5) cos^ + cos^)) 



U) 



The radiation pressure is due to the jump in photon linear momentum (pro- 
portional to the refractive index rij) when propagating from fluid 1 to fluid 2. 
T(8i,9t) is the Fresnel transmission coefficients of energy fluxes for circularly 



polarized beams (IChraibi et al\\2008a\ ) 



T(0i,0 t ) 



2(1 — 5) cos 9i cos 9t 



2(1 — 5) cos 6i cos i 



((1 - 5) cos 6i + cos t ) 2 (cos 6 t + (1 - 5) cos 6 t ) 2 



.2) 



where 6i and Qt are the incidence and transmission angles at the fluid interface 
respectively. They are such that 9i = arctan(^) and 9t = arcsin((l — 5) sin 0i) 
where 5 = (n 2 —fixj/n^ is the relative refractive index contrast. 
The fluid used in this experiment is a near-critical binary mixture in the two-phase 
regime that undergoes refractive index fluctuations of characteristic size £~, which 
increases as the critical point is neared. When the difference between the sam- 
ple temperature T and the critical temperature T c is large, £~ vanishes and each 
phas e at coexistence can b e considered as almost homogeneous (jCasner &: Delville! 
200ll : Casner et al. 20031 ) for the light-matter interaction. However, when T is 
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close to T c , £~ increases up to tens of nanometers at T-Tc=0.1K and the sam- 
ple significantly scatters light, producing strong linear momentum transfer to the 
flu id. Details on the ex p eriments will not be p resen ted h ere, however one can r efer 
to ISchroll et all (|2007l l. IWunenburger et al\ (|2010T ) and IChraibi et oil (|2010l ) for 



more details. 

The profiles are analytically predicted using the stream function formulation de- 
veloped in section 3 all along with the numerical solution of equation fj3. 16f> 
(considering h'(r = 2wo) = as observed in the experiments). 
Figure [17] presents comparisons between two-fluids numerical results, one-fluid 
analytical predictions and experimental steady interface profiles for two different 
experimental conditions (T-T c =l K, loq=5.3 fim, and T-T c =1.5 K, ujo=7.5 (j,m) 
and at four different values of the laser beam power. 

A very good agreement is obtained between the experimental data and the pre- 
dictions at each power value P = 154, 308, 462 and 616 mW. Insets of figure 
[T71 illustrate the linear power dependence of the deformation height due to the 
steady scattering flows disturbed by opposite radiation pressure effects. 
All these results validate both the numerical resolution and the theoretical model, 
which both represent important tools to describe these original physical phenom- 
ena. 
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Figure 17. Comparison of steady state interface profiles obtained from experiments, analytical 
solutions and numerical simulations. Top: T-T c =l K, ljo=5.3 /im, Bottom: T-T c =1.5 K, cjo=7.5 
fj,m. Insets show the variation of scattering deformation h max /u = h(r/u)Q = 2) as a function 
of beam power P (dot-dashed line: numerical, symbols: experiments). In all cases the optical 
wave propagates upward. Interface profiles were shifted vertically for clarity. 
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